rm(list=ls(all=TRUE))  # new
library(pscl)

load("rc.Rdata")
load("results.Rdata")
ideal.est = results[[1]]

### make a matrix with information on each legislator

trim <- function (x) gsub("^\\s+|\\s+$", "", x) # returns string w/o leading or trailing whitespace

make.members.matrix <- function(){
	members1 = colnames(ideal.est$x)
	members2 = strsplit(members1,split="[.]")
	for (i in 1:length(members2)){
		m = members2[[i]]
		mname = trim(paste(m[5:length(m)],collapse="."))
		members2[[i]] = c(m[1:4],mname)
	}
	t(matrix(unlist(members2),nrow=5))
}
members = make.members.matrix()
colnames(members)=c("type","cong","id","party","name")



### positions of members satisfying the criteria
get.pos <- function(type="",cong="",id="",party=""){
	cong1 = format(cong)
	id1 = format(id)
	pos1 = (members[,"type"] %in% type) | type==""
	pos2 = (members[,"cong"] %in% cong1) | cong==""
	pos3 = (members[,"id"] %in% id1) | id==""
	pos4 = (members[,"party"] %in% party) | party==""
	which(pos1&pos2&pos3&pos4)
}

### simulated ideal points of members satisfying the criteria
get.x <- function(type="",cong="",id="",party="",print.members=TRUE,save.att=FALSE){
	pos = get.pos(type,cong,id,party)
	att = members[pos,]
	if (print.members) {
		if (length(pos)==1) {
			print(cbind(pos,t(att)))
		} else { 
			print(cbind(pos,att))
		}
	}
	if (save.att){
		if (length(pos)==1){
			return(rbind(cong = att["cong"],type = att["type"],
			 party=att["party"],ideal = ideal.est$x[,pos,1]))

		} else {		
		return(rbind(cong = att[,"cong"],type = att[,"type"],
		 party=att[,"party"],ideal = ideal.est$x[,pos,1]))
		}
	} else{
		return(ideal.est$x[,pos,1])
	}
}


### search by name (partial match)
search.by.name <- function(mname){
	pos = grep(mname,members[,"name"])
	cbind(pos,members[pos,])
}


#### EXAMPLE ###################################################
#search.by.name("S[EU]T")  # search for SET or SUT
#search.by.name("S..TT")  # search for S??TT
#search.by.name("OSTETT")  # search for OSTETT
################################################################
#x1=get.x(id="29520")
#x105=get.x(id="29520",cong=105)
#x106=get.x(id="29520",cong=106)
#xHR103=get.x(cong=103,party="R",type="H")
#xP=get.x(type="P")
#x105_106 = get.x(id="29520", cong=c(105,106))
########################################################################

### posterior prob. that person "a" is more liberal than person "b" (A and B might 
###be the same person in different Congresses)

prob.more.lib <- function (a.id, a.cong, b.id, b.cong){
	a = get.x (id=a.id, cong=a.cong, print.members=FALSE)
	b = get.x (id=b.id, cong=b.cong, print.members=FALSE)
      r = list (dif.med = (median(b)-median(a)), prob = mean(a<b))
	r
}


### plot by id

plot.by.one.id <- function (id){ 
	pos = get.pos(id=id)
	cong = members [pos,"cong"]
	info.ideal = get.x(id = id )	
	colnames(info.ideal)=cong
	type = members [pos,"type"]
	

	info.ideal = info.ideal[,colnames(info.ideal)>101]
	info.ideal = info.ideal[,colnames(info.ideal)<112]

	boxplot(info.ideal, notch=TRUE, xlab="Congress",
    ylim=c(2.5*min(s.dem.med),2.5*max(s.rep.med)) )

	if (type[(length(pos))]=="S"){
		points(s.dem.med[c(-1,-12)],col=4,lwd=2,pch=23)#median of senate democrat
		points(s.rep.med[c(-1,-12)],col=2,lwd=2,pch=23)#median of senate republican
		legend("bottomleft",col=c(2,4),pch=23,pt.lwd=2,legend=c("Sen R","Sen D"), bty="n")
	} else {
		points(h.dem.med[c(-1,-12)],col=4,lwd=2,pch=23)#median of senate democrat
		points(h.rep.med[c(-1,-12)],col=2,lwd=2,pch=23)#median of senate republican
		legend("bottomleft",col=c(2,4),pch=23,pt.lwd=2,legend=c("House R","House D"), bty="n")

	}
}

plot.by.two.id <- function (id1,id2){ 
	
	pos1 = get.pos(id=id1)
	cong1 = members [pos1,"cong"]
	info.ideal1 = get.x(id = id1 )	

	pos2 = get.pos(id=id2)
	cong2 = members [pos2,"cong"]
	info.ideal2 = get.x(id = id2 )
	type = members [pos2,"type"]

	
	info.ideal = cbind(info.ideal1, info.ideal2)
	colnames(info.ideal)=c(cong1,cong2)
	
	info.ideal = info.ideal[,colnames(info.ideal)>101]
	info.ideal = info.ideal[,colnames(info.ideal)<112]
	boxplot(info.ideal, notch=TRUE, xlab="Congress",
    ylim=c(2.5*min(s.dem.med),2.5*max(s.rep.med)))

	if (type[(length(pos2))]=="S"){
		points(s.dem.med[c(-1,-12)],col=4,lwd=2,pch=23)#median of senate democrat
		points(s.rep.med[c(-1,-12)],col=2,lwd=2,pch=23)#median of senate republican
		legend("bottomleft",col=c(2,4),pch=23,pt.lwd=2,legend=c("Sen R","Sen D"), bty="n")
	} else {
		points(h.dem.med[c(-1,-12)],col=4,lwd=2,pch=23)#median of house democrat
		points(h.rep.med[c(-1,-12)],col=2,lwd=2,pch=23)#median of house republican
		legend("bottomleft",col=c(2,4),pch=23,pt.lwd=2,legend=c("House R","House D"), bty="n")
	}
	return(info.ideal)
}


# for table2 

ave.med <- function (id, type){ #take an average of medians for an individual	
	temp = get.x(id=id, type =type, print.members=FALSE)
	if (is.vector(temp)=="TRUE"){
		r = median(temp)
	} else {
		r.temp = rep(NA, times=ncol(temp))
		for (i in 1:ncol(temp)){
			r.temp[i]=median(temp[,i])
		}
		r = mean(r.temp)
	}
	r
}



compare.ave.med <- function (id){ # function for those who served in both chamber
	h.ave = ave.med (id, "H")
	s.ave = ave.med (id, "S")
	return = c(h.ave=h.ave, s.ave=s.ave)
}

##################################
# RUN
##################################

id.table = rownames( table(members[,"id"]) )

###################################################################
## change.mat is a matrix with posterior medians. 
## Each column from the third is for each congress, and each row is for each id. 
## For example, change.mat[1,3] shows the posterior median of 
#the first legislator in the 101st congress.
###################################################################
change.mat = matrix(NA, nrow = length(id.table), ncol = 14)
rownames(change.mat)=id.table
colnames(change.mat)=c("type","party",101:112)

for (i in 1:length(id.table)){
	temp = get.x(id=id.table[i],save.att=TRUE, print.members=FALSE)
	n.cong = dim(temp)[2]

	for(j in 1:n.cong){
		congr = as.numeric(temp["cong",j])-100
		change.mat[i,1]=temp["type",j]
		change.mat[i,2]=temp["party",j]
		change.mat[i,(congr+2)] = median(as.numeric(temp[c(-1,-2,-3),j]),na.rm = TRUE)	
	}
}	

change.mat.median = change.mat[,c(-1,-2)]


################
#   FIGURE 5   #
################

par(mfrow=c(1,2))

####### (1) HOUSE
h.dem.vec = (change.mat[,"type"]=="H" & change.mat[,"party"]=="D")
h.rep.vec = (change.mat[,"type"]=="H" & change.mat[,"party"]=="R")
matplot(102:111,t(change.mat.median[h.dem.vec,c(-1,-1*ncol(change.mat.median))]),type='l',ylim=c(-3,3),
col=8, main="House", ylab= expression(paste(plain('"Anchors Away"'),plain('ideal point estimate'))), xlab="Congress")
matlines(102:111,t(change.mat.median[h.rep.vec,c(-1,-1*ncol(change.mat.median)) ]),type='l',ylim=c(-3,3),col=1)


####### (2) SENATE
s.dem.vec = (change.mat[,"type"]=="S" & change.mat[,"party"]=="D")
s.rep.vec = (change.mat[,"type"]=="S" & change.mat[,"party"]=="R")

matplot(102:111,t(change.mat.median[s.dem.vec,c(-1,-1*ncol(change.mat.median)) ]),type='l',
ylim=c(-3,3),col=8, main="Senate", ylab= expression(paste(plain('"Anchors Away"'),plain('ideal point estimate'))), xlab="Congress")
matlines(102:111,t(change.mat.median[s.rep.vec,c(-1,-1*ncol(change.mat.median)) ]),type='l',ylim=c(-3,3),col=1)

legend("topleft", col=1, legend=c("Rep" ), box.col=0, lty=1, lwd=2)
legend("topright", col=8, legend=c("Dem" ), box.col=0, lty=1, lwd=2)


########################################
# Table 2 Senators who served in the House    #
########################################

p6.ids = c(14400,15116,14809,14803,15407,15011,14812,14661,14871,
29367,15429,15424,15020,29141,15633,29523,29108,15425,15021,29142,29306,
15071,14852,14858,29345,15406,14651,29537,29148,15015,39310,29732,29512,
29369,29740,29566,29754,29936,29918,29373,29555,29548,29909,15408,29147,
29389,29534,29924,29906,20115,20730,29935,29735,20101,29722,29386)

p6.names = rep(NA, length(p6.ids))
compare.HS = matrix(NA, nrow = length(p6.ids), ncol=2)

#memo: 95407 campbell (excluded after changing party affiliation)

for ( i in 1:length(p6.ids)){
	p6.names[i] = members[get.pos(id=p6.ids[i])[1],"name"]
	compare.HS [i,1] = compare.ave.med(p6.ids[i])["h.ave"]
	compare.HS [i,2] = compare.ave.med(p6.ids[i])["s.ave"]
}
p6.table = cbind(p6.names,round(compare.HS,digit=3))
rownames(p6.table)=p6.ids
colnames (p6.table)=c("NAME","H.AVE","S.AVE")
p6.table
xs=read.csv("controls_t2.csv")
pos1 = match(p6.ids[1],xs[,1])
xs.new = xs[pos1,]

for (i in 2:length(p6.ids)){
print(i)

	pos = match(p6.ids[i],xs[,1])
	xs.new=rbind(xs.new, xs[pos,])
}

colnames(xs.new)=colnames(xs)


HSdifference = as.numeric(p6.table[,3])-as.numeric(p6.table[,2])  
t2 = lm(HSdifference~xs.new[,2]+xs.new[,3]+xs.new[,4]+xs.new[,5]+xs.new[,6]+xs.new[,7])
summary(t2)

#########################################
# footnote 12: brigde decision sponsors vs. non-bridge bill sponsors
#########################################
library(foreign)
id.cosponsor = read.dta("sponsorICPSR121.dta")

id.list.nonbridge = id.cosponsor [(id.cosponsor[, "cosponsorshipdecisionvote" ]==0),c(  "sponsoricpsr" ,"congress")]
id.list.bridge = id.cosponsor [(id.cosponsor[, "cosponsorshipdecisionvote" ]==1), c(  "sponsoricpsr" ,"congress")]

nonbridge.ideal = rep(NA, times=dim(id.list.nonbridge)[1])
bridge.ideal = rep(NA, times=dim(id.list.bridge)[1])

for (i in 1:length(nonbridge.ideal)){
	temp=	(rownames(change.mat)==id.list.nonbridge[i, "sponsoricpsr"])
	if (sum(temp,na.rm=TRUE)==0){
		nonbridge.ideal[i]=c(NA)
	} else {
		nonbridge.ideal[i]=change.mat[temp,(id.list.nonbridge[i,"congress"]-98)]
	}
}
for (i in 1:length(bridge.ideal)){
	temp=	(rownames(change.mat)==id.list.bridge[i, "sponsoricpsr"])
	if (sum(temp,na.rm=TRUE)==0){
		bridge.ideal[i]=c(NA)
	} else {
		bridge.ideal[i]=change.mat[temp,(id.list.bridge[i,"congress"]-98)]
	}
}

par(mfrow=c(1,2))
hist(as.numeric(nonbridge.ideal),freq = FALSE,xlim=c(-3,3),main="Non-bridge bill sponsors",
xlab=expression(paste(plain('"Anchors Away"'),plain('ideal point estimate'))))
hist(as.numeric(bridge.ideal),freq = FALSE,xlim=c(-3,3), main="bridge bill sponsors",
xlab=expression(paste(plain('"Anchors Away"'),plain('ideal point estimate'))))

mean(as.numeric(bridge.ideal),na.rm=TRUE)
mean(as.numeric(nonbridge.ideal),na.rm=TRUE)
sd(as.numeric(bridge.ideal),na.rm=TRUE)
sd(as.numeric(nonbridge.ideal),na.rm=TRUE)
t.test(x=as.numeric(nonbridge.ideal),y=as.numeric(bridge.ideal))


##########
# footnote 12: partisanship scores
###########

pv.bridge = id.cosponsor [(id.cosponsor[, "cosponsorshipdecisionvote" ]==1&
	id.cosponsor[, "unanimous" ]==0),"partisanshipofvote"]
pv.nonbridge = id.cosponsor [(id.cosponsor[, "cosponsorshipdecisionvote" ]==0&
	id.cosponsor[, "unanimous" ]==0),"partisanshipofvote"]

mean(pv.bridge)
mean(pv.nonbridge)

t.test(x=pv.bridge, y=pv.nonbridge)

par(mfrow=c(1,2))
hist(pv.bridge, freq=FALSE, main="Partisanship scores for bridge bills",xlab="Partisanship score")
hist(pv.nonbridge,freq=FALSE, main="Partisanship scores for non-bridge bills",xlab="Partisanship score")



##############
# Table 3: party switchers
#################

# posterior for median
h.dem.meds = matrix (NA, ncol=12, nrow = dim(ideal.est$x)[1])
h.rep.meds = matrix (NA, ncol=12, nrow = dim(ideal.est$x)[1])
s.dem.meds = matrix (NA, ncol=12, nrow = dim(ideal.est$x)[1])
s.rep.meds = matrix (NA, ncol=12, nrow = dim(ideal.est$x)[1])
colnames(h.dem.meds)=c(101:112)
colnames(h.rep.meds)=c(101:112)
colnames(s.dem.meds)=c(101:112)
colnames(s.rep.meds)=c(101:112)


for (i in 1:12){

	temp1 = get.x (type="H", cong = (i+100), party ="D", print.members=FALSE)
	temp2 = get.x (type="H", cong = (i+100), party ="R", print.members=FALSE)
	temp3 = get.x (type="S", cong = (i+100), party ="D", print.members=FALSE)
	temp4 = get.x (type="S", cong = (i+100), party ="R", print.members=FALSE)


	for (j in 1:nrow(h.dem.meds)){

		h.dem.meds [j,i] = median(as.numeric(temp1[j,])) 
		h.rep.meds[j,i] = median(as.numeric(temp2[j,]))
		s.dem.meds[j,i] = median(as.numeric(temp3[j,]))
		s.rep.meds[j,i] = median(as.numeric(temp4[j,]))
	}
	
}
#party median
h.dem.med = rep(NA,times=12)
h.rep.med = rep(NA,times=12)
s.dem.med = rep(NA,times=12)
s.rep.med = rep(NA,times=12)

for (i in 1:length(h.dem.med)){
	h.dem.med[i]= median(h.dem.meds[,i])
	h.rep.med[i]= median(h.rep.meds[,i])
	s.dem.med[i]= median(s.dem.meds[,i])
	s.rep.med[i]= median(s.rep.meds[,i])
}

#version1(difference in the distance from a party median
dif.dis.from.pmed <- function (ideal.mat, med.mat){

	dis = matrix(NA, nrow=nrow(ideal.mat), ncol = ncol(ideal.mat))
	colnames(dis)=colnames(ideal.mat)
	for (i in 1:ncol(ideal.mat)){
		dis[,i] = med.mat[,(as.numeric(colnames(ideal.mat)[i])-100)]-ideal.mat[i]
	}

	dif.dis.mat = matrix(NA, nrow=nrow(ideal.mat), ncol = (ncol(ideal.mat)-1))
	colnames(dif.dis.mat)=colnames(ideal.mat)[-1]
	dif.dis.med = rep(NA, times = ncol(dif.dis.mat))
	names(dif.dis.med)=colnames(dif.dis.mat)
	dif.prob = rep(NA, times = ncol(dif.dis.mat))

	for (i in 1:ncol(dif.dis.mat)){
		dif.dis.mat[,i]= dis[,(i+1)]-dis[,i]
		dif.dis.med[i] = median(dif.dis.mat[,i])
		dif.prob[i] = mean(dif.dis.mat[,i]>0)
	}
	r = list(dif.dis.mat=dif.dis.mat, dif.dis.med =dif.dis.med, dif.prob = dif.prob)
	r
} 

#version2
dif.id <- function (ideal.mat){

	dif.mat = matrix(NA, nrow=nrow(ideal.mat), ncol = (ncol(ideal.mat)-1))
	colnames(dif.mat)=colnames(ideal.mat)[-1]
	dif.med = rep(NA, times = ncol(dif.mat))
	names(dif.med)=colnames(dif.mat)  #ith cong-(i-1)th cong
	dif.prob = rep(NA, times = ncol(dif.mat))

	for (i in 1:ncol(dif.mat)){
		dif.mat[,i]= ideal.mat[,(i+1)]-ideal.mat[,i]
		dif.med[i] = median(dif.mat[,i])
		dif.prob[i] = mean(dif.mat[,i]>0)
	}
	r = list(dif.mat=dif.mat, dif.med =dif.med, dif.prob = dif.prob)
	r
}



#### 
#a) SHELBY
shelby = plot.by.two.id(14659,94659)
#shelby.dif = dif.dis.from.pmed (shelby,s.dem.meds)
shelby.dif.id = dif.id(shelby)
shelby.dif.id$dif.med
shelby.dif.id$dif.prob

#b)CAMPBELL
campbell = plot.by.two.id(15407,95407)[,-1] #first col is for House
campbell.dif = dif.dis.from.pmed (campbell,s.dem.meds)#[-1] 

campbell.dif.id = dif.id(campbell)
round (campbell.dif.id$dif.med,digit=2)
round(campbell.dif.id$dif.prob,digit=2)

#c) DEAL (gradually adjusted?)
deal = plot.by.two.id(29342,99342)
deal.dif.id = dif.id(deal)
deal.dif.id[-1]
round(deal.dif.id[[2]],digit=2)

#D) h.HAYES
hayes = plot.by.two.id(15418,95418)
hayes.dif.id = dif.id(hayes)
hayes.dif.id[-1]
round(hayes.dif.id[[2]],digit=2)

#E) LAUGHLIN GREG
laughlin = plot.by.two.id(15611,95611)

laughlin.dif.id = dif.id(laughlin)
laughlin.dif.id[-1]
round(laughlin.dif.id[[2]],digit=2)


#F) M.Parker

parker = plot.by.two.id(15617,95617)

parker.dif.id = dif.id(parker)
parker.dif.id[-1]
round(parker.dif.id[[2]],digit=2)


#G) TAUZIN BILLY
tauzin = plot.by.two.id(14679,94679)

tauzin.dif.id = dif.id(tauzin)
tauzin.dif.id[-1]
round(tauzin.dif.id[[2]],digit=2)


#H)FORBES 
forbes = plot.by.two.id(29542,99542)

forbes.dif.id = dif.id(forbes)
forbes.dif.id[-1]
round(forbes.dif.id[[2]],digit=2)


#I)  GOODE (D to I)
goode1 = plot.by.two.id(29767,99767)

goode1.dif.id = dif.id(goode1)
goode1.dif.id[-1]
round(goode1.dif.id[[2]],digit=2)


#J)MARTINEZ MATTHEW G
martinez = plot.by.two.id(14879,94879)

martinez.dif.id = dif.id(martinez)
martinez.dif.id[-1]
round(martinez.dif.id[[2]],digit=2)


#K) GOODE (I to R)
goode2 =  plot.by.two.id(99767,89767)

goode2.dif.id = dif.id(goode2)
goode2.dif.id[-1]
round(goode2.dif.id[[2]],digit=2)


#L) JEFFORDS
jeffords = plot.by.two.id (14240,94240)

jeffords.dif.id = dif.id(jeffords)
jeffords.dif.id[-1]
round(jeffords.dif.id[[2]],digit=2)


#M) HALL 
hall= plot.by.two.id(14828,94828)

hall.dif.id = dif.id(hall)
hall.dif.id[-1]
round(hall.dif.id[[2]],digit=2)


#N) ALEXANDER
alexander = plot.by.two.id(20327,90327)

alexander.dif.id = dif.id(alexander)
alexander.dif.id[-1]
round(alexander.dif.id[[2]],digit=2)

#O)GRIFFITH
griffith = plot.by.two.id(20901,90901)

griffith.dif.id = dif.id(griffith)
griffith.dif.id[-1]
round(griffith.dif.id[[2]],digit=2)

#P)SPECTER
specter = plot.by.two.id(14910,94910)

specter.dif.id = dif.id(specter)
specter.dif.id[-1]
round(specter.dif.id[[2]],digit=2)
